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Abstract. An introduction to the Propp-Wilson method of coupling- from-the- 
past for the Ising model is presented. It enables one to obtain exact samples from 
the equilibrium spin distribution for ferromagnetic interactions. Both uniform and 
random quenched magnetic fields are included. The time to couple from the past 
and its standard deviation are shown as functions of system size, temperature, and 
the field strength. The time should be an estimate of the ergodic time for the Ising 
system. 

The standard importance-sampling Monte Carlo algorithm jj| for statics 
has a number of difficulties. One is its inability to ensure that the generated 
spin configurations are uncorrelated and are representative of the equilib- 
rium distribution. The standard way to overcome this difficulty is to throw 
away a large number of Monte Carlo steps between analyzed configurations. 
Of course this decreases the number of configurations obtained in a given 
amount of computer time. In addition, there can remain the question of 
whether the generated spin configurations are indeed drawn from the un- 
derlying equilibrium distribution. For example, consider a spin-glass model 
at low temperatures, where ergodicity has been broken for reasonable sim- 
ulation times and the system might have many ground states and many 
more metastable states. Then it is extremely difficult to decide whether the 
standard Monte Carlo method obtains spin configurations drawn from the 
underlying equilibrium spin distribution. 

In this short paper the square-lattice Ising model with nearest-neighbor 
ferromagnetic (J>0) interactions and periodic boundary conditions is con- 
sidered. The Hamiltonian is TL = — J^luj) (J i (J j ~ S» Hi&i, with cr^=±l, the 
first sum over nearest-neighbor spins, and Hi is the local field. Two cases are 
considered, uniform field Hi=H, and quenched random field with Hi cho- 
sen to be ±iJ with equal probability. We take J=l, which makes the exact 
critical temperature T c ?»2.26 • • •. 

One Monte Carlo step consists of two procedures. First randomly choose 
a spin i from the LxL lattice. The probability of Gi being set to +1 in the 
next time step is 



Prob(cr i = l) = < 1 + cxp 



-2J ^2 °j ~ 2H ' 



]=nn 
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where the first sum is over the nearest-neighbors of spin i and T is the tem- 
perature. The second step is to use a uniformly distributed random number r, 
and set <7i=+l in the next time step if r<Prob(<7i = 1); otherwise set a——!. 
This is the single-spin flip heat-bath dynamic, which for the Ising model is 
the Glauber dynamic. I use the Ziff four-tap random number generator ||. 
It is important to note that we do not consider the probability that a spin 
flips, but rather the probability of the spin being up, i.e. +1. 

Recently researchers in the statistics community have devised methods 
that can guarantee that generated spin configurations from a Monte Carlo 
simulation are drawn from the underlying equilibrium distribution. There the 
study goes under the name of Markov random fields || , which should not be 
confused with the random fields in the Hamiltonian. The method used here is 
the 1996 method of Propp and Wilson It is referred to as ' coupling- from- 
the-past' or the 'Propp- Wilson method'. Here we bring the single-spin-flip 
version of the Propp- Wilson method to the attention of practitioners in the 
statistical mechanical community. We obtain a time for coupling from the 
past, which we call the Propp-Wilson time, tpw We will show how ipw 
changes with T, L, H, and the randomness in H. The time tpw can be 
viewed as an upper bound to the ergodic time r c of the simulation [Q. In 
the spirit of a workshop, this study was only initiated about one month prior 
to the workshop. Consequently, it should be viewed as an exploratory and 
preliminary work. 

The basic idea of the Propp-Wilson method is to not run the Markov 
chain forward in time (which all normal importance-sampling Monte Carlo 
procedures do), but rather to run the Markov chain from a time —t in the 
past up to t—0. Typically it is not known how to go from a spin configuration 
at one time to one at an earlier time. However, given a spin configuration at 
time — t, it is easy to propagate it to time t=0 — just use the algorithm for 
forward propagation. What is the difference between going from —t to or 
going from to tl It lies only in our knowledge of the sequence of random 
numbers used. At t=0 one has no idea which spin will be picked next or which 
random number the spin-up probability will be compared with. However, at 
all times in the past, both the spin picked and the random number to compare 
with are known. 

Now consider starting two spin configurations at time — t, and propagating 
forward in time toward t=0. The same random numbers are used for each 
spin configuration, since the same history of random numbers applies for 
times t<0. At some time —t c these two spin configurations might coalesce, 
i.e. they both become the same spin configuration at time —t c . From this 
time onward to time t=0 they will remain coalesced, since the same random 
numbers propagate them forward to t=0. This coalescence gives the Propp- 
Wilson method the 'coupling from the past' name. 

Here is the Propp-Wilson prescription for simulation of N Ising spins. 
First obtain a very long history of randomly picked spins and random num- 
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bers used to compare against Prob(<7i = +1). Start with all 2 N possible spin 
configurations at time t=—l. Propagate all configurations forward in time 
to t=0. Remember, this is done using the same set of random numbers for 
each configuration. If all 2 N configurations have coalesced to the same spin 
configuration at t=0, this is the spin configuration you want. If not, start all 
2 N spin configurations at t=— 2, propagate forward with the same random 
number history, and see if they have coalesced at i=0. This is continued until 
at some time — tpw the coalescence of all 2 N spin configurations has occurred. 
The spin configuration obtained in this fashion is guaranteed to be an exact 
sampling of the equilibrium spin distribution. An exact sample means that 
the spin configuration is drawn according to the probability associated with 
the equilibrium distribution. Note that care must be taken to ensure that the 
long random-number history is never exceeded before the coalescence occurs. 

Why does this give an exact sampling? Let us imagine that we have run 
an infinite ensemble of spin configurations from t—— oo to t=— tpw using 
different random numbers for each. Since we have run our Markov chain 
infinitely long, we are guaranteed that the spin configurations at t=— ipw 
are distributed according to the equilibrium spin distribution. Each of the 
starting spin configurations will have evolved into one of the 2 N spin config- 
urations at t=— tpw- No matter which spin configuration the system is in at 
time — tpw, we know that at time t=0 all 2 N configurations have coalesced to 
a single spin configuration since we have used our history of random numbers 
which is the same for each configuration to run from t=— tpw to t=0. Thus 
the sample generated is drawn from the equilibrium spin distribution. 

In principle only a coupling-from-the-past method, like the Propp- Wilson 
method, is required to obtain a perfectly random sample. In practice another 
relationship is required since the number of states grows very quickly with 
lattice size. If all 2 N states could be handled, the partition function could 
be calculated exactly. For special cases it is possible to introduce a partial 
ordering on the states. For the Ising model, we say that two spin configu- 
rations x and y have a partial ordering, x ^ y, if for each lattice site the 
spin in x is down whenever the spin in y is down. The transition probability 
of (Q) preserves the partial ordering for ferromagnetic interactions. Conse- 
quently, rather than insist that all 2 N configurations have coalesced, only the 
two extreme partial ordering states, all spins up and all spins down, need to 
be considered. If these two states have coalesced, then all other states have 
coalesced. It is this desired property of partial ordering that restricts one to 
ferromagnetic interactions. 

The Propp- Wilson method was run to obtain tpw as a function of L, T, 
and H . For each value of parameters 1000 exact samples were obtained. Note 
that if coalescence occurs at time — ipw> it also occurs for any earlier time. 
This allows a bisection method to be used to find tpw- Figure 1(a) shows 
the average value of the Propp- Wilson time, [tpw], m Monte Carlo Steps per 
Spin (MCSS), at H=0. Since the phase space grows with system size, [tpw] 
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Fig. 1. (a) The average Propp- Wilson time, [tpw], is shown as a function of tem- 
perature for H=0 and L=4, 8, and 16. The vertical line is the exact value of the 
critical temperature, T c . The lines are guides for the eye. (b) The L dependence 
on a log-log plot is shown for [tpw] (filled symbols joined by solid lines) and the 
standard deviation of tpw (corresponding open symbols joined by dashed lines). 
Both of these grow faster than power laws below T c , slower than power laws above 
T c , and are expected to grow as a power law at T c . 



increases with L. Furthermore, since it takes longer to sample all of phase 
space at lower temperatures, [ip\y] increases as T is lowered. 

For a d dimensional system with correlation length £, below but near 
T c the ergodic time is expected to scale as r c ocL z exp [const (L/^) d_1 ] 
Here z is the exponent associated with critical slowing down. For T>T C the 
'intrinsic' order parameter relaxation time r is the longest relaxation time 
in the system Jjl], so above T c one expects toc£ z cx|1 — T/T c \~ vz where v is 
the correlation length exponent. Right at T c this yields tozL z . Figure 1(b) 
shows the L dependence of [ipw] for H—0 and various temperatures. The 
graph illustrates that this time grows slower (faster) than a power law above 
(below) T c . At T c this time seems to scale as a power law. However, care 
must be used in associating [ipw] with either t or r . One might expect that 
[^pw] yields an upper bound on the larger of r or r. As seen in Fig. 1 even at 
high temperatures [ipw] increases with L. These additional size dependent 
factors in [tpw] might also be expected to be present at lower temperatures. 
In a recent thesis || [ipw]oci/ 2 was postulated, and the value X?«2.44 was 
found with system sizes L up to 310. (This assumes that the presented values 
are at T c , which was not stated explicitly.) This is somewhat high compared 
to the accepted value of z«2.21 [||. It is also possible to calculate higher 
moments of tpw- Figure 1(b) shows the standard deviation per spin for tpw 
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Fig. 2. The average Propp-Wilson time, [tpw], in MCSS is shown at T—2.0 as 
a function of H. The lines are guides for the eye. Note that the points at H=0 
are the same for both plots, (a) Is for a uniform magnetic field. For L—16, H=0.1 
is the smallest value shown, (b) Is for quenched randomly distributed fields equal 
to —H or +H. Four different distributions of fields are shown for L=4, and one 
distribution for L=8 and L—16. For L=16, H=1.5 is the smallest value shown. 



as a function of L. The behavior seems to be close to those of [ipw]- Note 
that the standard deviation becomes very large, which necessitates using a 
large number of samples to obtain a reasonable value for [ip\y]- 

Figure 2(a) shows [tpw] a t T—2.0 for different values of L and uniform 
applied field, H. At strong fields the larger system sizes have larger values 
of [f pw] i for intermediate fields all system sizes have similar values of [ipw] , 
and for weak fields the smaller system sizes break away from this curve and 
go to their maximum at H=0. This behavior is reminiscent of the cross-overs 
encountered in the study of metastability, where for T<T C as the field is 
decreased different regimes are encountered in a finite system (7). 

To illustrate the application of the Propp-Wilson method to a random- 
field Ising model, Fig. 2(b) shows data for [tpw] as a function of the strength 
of the random field for various quenched field configurations for L=A, and one 
quenched field configuration for L—8 and 16. The trends seem to be similar 
to those of Fig. 2(a). However the values for [£p\y] seem to increase faster 
as H decreases than in the uniform field case. This is in agreement with 
the physical idea that a system with random quenched disorder has a longer 
ergodicity time than one without disorder. For more discussion of random 
fields see references [isLEillTT 
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In summary, it is possible to utilize coupling-from-the-past ideas to ob- 
tain perfect samples of Ising spin configurations from the equilibrium spin 
distribution. To use the ideas in practice, however, another relationship is 
essential, that of a partial ordering preserved by the transition probabilities. 
This restriction allows one to study random fields, but no partial ordering is 
known for systems with antiferromagnetic bonds. The value for the average 
and standard deviation for tpw seem to have reasonable dependences on sys- 
tem size, field strength, and temperature if ip\y is viewed as an upper bound 
for the ergodic time of the system. 

It is possible to use cluster dynamics instead of single-spin- flip dynamics. 
Since this decreases the ergodic time near T c , this allows one to obtain per- 
fect samples for large systems — so far the largest single example I have seen 
had L=2100 in a preprint by Propp and Wilson. It has also been shown that 
the expected running time for the Propp- Wilson method should be within a 
constant multiple of the maximum expected time for two states to coalesce in 
the Markov chain, and that this running time is near optimal for exact sam- 
pling. Thus the ability to obtain results for large systems at temperatures 
above and near T c for models that have a partial ordering seems realistic. 
Unfortunately, these partial orderings are extremely difficult to find in more 
general models. Also, the long expected run times seem to make the applica- 
tion of these types of algorithms rather limited in studies of low-temperature 
behavior for systems with quenched randomness and high barriers. 

This work is supported in part by the NSF through grant number 9871455, 
and through the Supercomputer Computations Research Institute which is 
funded by the U.S. DoE and the State of Florida. 
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